Load the necessary libraries
library(car) #for regression diagnostics
library(broom) #for tidy output
library(ggfortify) #for model diagnostics
library(sjPlot) #for outputs
library(knitr) #for kable
library(effects) #for partial effects plots
library(emmeans) #for estimating marginal means
library(ggeffects) #for plotting marginal means
library(MASS) #for glm.nb
library(MuMIn) #for AICc
library(tidyverse) #for data wrangling
library(modelr) #for auxillary modelling functions
library(DHARMa) #for residual diagnostics plots
library(performance) #for residuals diagnostics
library(see) #for plotting residuals
library(patchwork) #grid of plots
library(scales) #for more scales
An ecologist studying a rocky shore at Phillip Island, in southeastern Australia, was interested in how clumps of intertidal mussels are maintained (Quinn 1988). In particular, he wanted to know how densities of adult mussels affected recruitment of young individuals from the plankton. As with most marine invertebrates, recruitment is highly patchy in time, so he expected to find seasonal variation, and the interaction between season and density - whether effects of adult mussel density vary across seasons - was the aspect of most interest.
The data were collected from four seasons, and with two densities of adult mussels. The experiment consisted of clumps of adult mussels attached to the rocks. These clumps were then brought back to the laboratory, and the number of baby mussels recorded. There were 3-6 replicate clumps for each density and season combination.
Format of quinn.csv data files
| SEASON | DENSITY | RECRUITS | SQRTRECRUITS | GROUP |
|---|---|---|---|---|
| Spring | Low | 15 | 3.87 | SpringLow |
| .. | .. | .. | .. | .. |
| Spring | High | 11 | 3.32 | SpringHigh |
| .. | .. | .. | .. | .. |
| Summer | Low | 21 | 4.58 | SummerLow |
| .. | .. | .. | .. | .. |
| Summer | High | 34 | 5.83 | SummerHigh |
| .. | .. | .. | .. | .. |
| Autumn | Low | 14 | 3.74 | AutumnLow |
| .. | .. | .. | .. | .. |
| SEASON | Categorical listing of Season in which mussel clumps were collected independent variable |
| DENSITY | Categorical listing of the density of mussels within mussel clump independent variable |
| RECRUITS | The number of mussel recruits response variable |
| SQRTRECRUITS | Square root transformation of RECRUITS - needed to meet the test assumptions |
| GROUPS | Categorical listing of Season/Density combinations - used for checking ANOVA assumptions |
Mussel
quinn = read_csv('../public/data/quinn.csv', trim_ws=TRUE)
glimpse(quinn)
## Rows: 42
## Columns: 5
## $ SEASON <chr> "Spring", "Spring", "Spring", "Spring", "Spring", "Sprin…
## $ DENSITY <chr> "Low", "Low", "Low", "Low", "Low", "High", "High", "High…
## $ RECRUITS <dbl> 15, 10, 13, 13, 5, 11, 10, 15, 10, 13, 1, 21, 31, 21, 18…
## $ SQRTRECRUITS <dbl> 3.872983, 3.162278, 3.605551, 3.605551, 2.236068, 3.3166…
## $ GROUP <chr> "SpringLow", "SpringLow", "SpringLow", "SpringLow", "Spr…
summary(quinn)
## SEASON DENSITY RECRUITS SQRTRECRUITS
## Length:42 Length:42 Min. : 0.00 Min. :0.000
## Class :character Class :character 1st Qu.: 9.25 1st Qu.:3.041
## Mode :character Mode :character Median :13.50 Median :3.674
## Mean :18.33 Mean :3.871
## 3rd Qu.:21.75 3rd Qu.:4.663
## Max. :69.00 Max. :8.307
## GROUP
## Length:42
## Class :character
## Mode :character
##
##
##
Since we intend to model both SEASON and DENSITY as categorical variables, we need to explicitly declair them as factors.
quinn = quinn %>% mutate(SEASON = factor(SEASON, levels=c('Spring', 'Summer', 'Autumn', 'Winter')),
DENSITY = factor(DENSITY))
Model formula: \[ \begin{align} y_i &\sim{} \mathcal{Pois}(\lambda_i)\\ ln(\mu_i) &= \boldsymbol{\beta} \bf{X_i}\\[1em] \end{align} \]
where \(\boldsymbol{\beta}\) is a vector of effects parameters and \(\bf{X}\) is a model matrix representing the intercept and effects of season, density and their interaction on mussel recruitment.
ggplot(quinn, aes(y=RECRUITS, x=SEASON, fill=DENSITY)) +
geom_boxplot()
Conclusions:
Lets mimic the effect of using a log link, by using log scaled y-axis.
ggplot(quinn, aes(y=RECRUITS, x=SEASON, fill=DENSITY)) +
geom_boxplot() +
scale_y_log10()
Conclusions:
We actually have a couple of initial options:
log transform the response variable (RECRUITS) and fit against a Gaussian distribution. If we do so, we would need to ammend the response for cases when the response is zero (as the log of 0 is not defined). We could simply add 1 to each count before log transformed.
fit against a Poisson distribution with a log link
I will include this, just for an illustration on how to fit such a model, we will not persue it further.
quinn.glmG <- glm(log(RECRUITS+1) ~ DENSITY*SEASON, data=quinn, family=gaussian)
quinn.glm <- glm(RECRUITS ~ DENSITY*SEASON, data=quinn,
family=poisson(link='log'))
autoplot(quinn.glm,which=1:6)
Conclusions:
performance::check_model(quinn.glm)
performance::check_overdispersion(quinn.glm)
## # Overdispersion test
##
## dispersion ratio = 3.309
## Pearson's Chi-Squared = 112.497
## p-value = < 0.001
performance::check_zeroinflation(quinn.glm)
## # Check for zero-inflation
##
## Observed zeros: 2
## Predicted zeros: 0
## Ratio: 0.00
Conclusions:
quinn.resid = simulateResiduals(quinn.glm, plot=TRUE)
testResiduals(quinn.resid)
## $uniformity
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.2008, p-value = 0.06761
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 1.7813, p-value < 2.2e-16
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 6, observations = 42, p-value < 2.2e-16
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.00000000 0.04761905
## sample estimates:
## outlier frequency (expected: 0.00880952380952381 )
## 0.1428571
## $uniformity
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.2008, p-value = 0.06761
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 1.7813, p-value < 2.2e-16
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 6, observations = 42, p-value < 2.2e-16
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.00000000 0.04761905
## sample estimates:
## outlier frequency (expected: 0.00880952380952381 )
## 0.1428571
testDispersion(quinn.resid)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 1.7813, p-value < 2.2e-16
## alternative hypothesis: two.sided
testZeroInflation(quinn.resid)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 10.87, p-value = 0.016
## alternative hypothesis: two.sided
#testTemporalAutocorrelation(quinn.glm1)
Conclusions:
## goodness of fit
1-pchisq(quinn.glm$deviance, df=quinn.glm$df.residual)
## [1] 2.457479e-12
## any evidence of overdispersion
quinn.glm$deviance/quinn.glm$df.residual
## [1] 3.677366
Conclusions:
quinn %>% group_by(SEASON, DENSITY) %>%
summarise(Zeros=sum(RECRUITS==0),
Prop=Zeros/n(),
Mean=mean(RECRUITS))
x=rpois(100000,lambda=2.67)
tab.1 = table(x==0)
tab.1/sum(tab.1)
##
## FALSE TRUE
## 0.93188 0.06812
## OR, over the entire data
## is this due to excessive zeros (zero-inflation)
tab=table(quinn$RECRUITS==0)
tab/sum(tab)
##
## FALSE TRUE
## 0.95238095 0.04761905
## 5% is not many.. so it cant be zero-inflated
## how many 0's would we expect from a poisson distribution with a mean similar to our mean
mean(quinn$RECRUITS)
## [1] 18.33333
x=rpois(100000,lambda=mean(quinn$RECRUITS))
tab.1 = table(x==0)
tab.1/sum(tab.1)
##
## FALSE
## 1
Conclusion:
In light of the discovery that the Poisson model was overdispersed, we will explore a negative binomial model. Rather than assume that the variance is equal to the mean (dispersion of 1), the dispersion parameter is estimated. Negative binomial models are also able to account for some level of zero-inflation.
library(MASS)
quinn.nb = glm.nb(RECRUITS ~ DENSITY*SEASON, data=quinn)
autoplot(quinn.nb,which=1:6)
performance::check_model(quinn.nb)
quinn.resid = simulateResiduals(quinn.nb, plot=TRUE)
testResiduals(quinn.resid)
## $uniformity
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.11213, p-value = 0.626
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 0.89834, p-value = 0.624
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 1, observations = 42, p-value = 0.66
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.00000000 0.04761905
## sample estimates:
## outlier frequency (expected: 0.0102380952380952 )
## 0.02380952
## $uniformity
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.11213, p-value = 0.626
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 0.89834, p-value = 0.624
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 1, observations = 42, p-value = 0.66
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.00000000 0.04761905
## sample estimates:
## outlier frequency (expected: 0.0102380952380952 )
## 0.02380952
testDispersion(quinn.resid)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 0.89834, p-value = 0.624
## alternative hypothesis: two.sided
testZeroInflation(quinn.resid)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 6.1728, p-value = 0.104
## alternative hypothesis: two.sided
## goodness of fit
1-pchisq(quinn.nb$deviance, df=quinn.nb$df.residual)
## [1] 0.01311451
## any evidence of overdispersion
quinn.nb$deviance/quinn.nb$df.residual
## [1] 1.614218
AICc(quinn.glm, quinn.nb)
Conclusions:
Conclusions:
plot_model(quinn.nb, type='eff', terms=c('SEASON', 'DENSITY'))
plot(allEffects(quinn.nb),multiline=TRUE, ci.style='bar')
plot(allEffects(quinn.nb),multiline=TRUE, ci.style='bar', type='link')
ggpredict(quinn.nb, c('SEASON', 'DENSITY')) %>% plot()
ggemmeans(quinn.nb, ~SEASON*DENSITY) %>% plot()
Lets start with the estimated coeffcicients on the link (log) scale
summary(quinn.nb)
##
## Call:
## glm.nb(formula = RECRUITS ~ DENSITY * SEASON, data = quinn, init.theta = 9.022467857,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.8704 -0.8274 0.0000 0.4999 2.1866
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.3026 0.1875 12.283 < 2e-16 ***
## DENSITYLow 0.1133 0.2742 0.413 0.67934
## SEASONSummer 1.5721 0.2389 6.581 4.69e-11 ***
## SEASONAutumn 0.6763 0.2492 2.714 0.00664 **
## SEASONWinter -0.5680 0.2881 -1.971 0.04870 *
## DENSITYLow:SEASONSummer -0.8970 0.3509 -2.556 0.01059 *
## DENSITYLow:SEASONAutumn -0.1881 0.3788 -0.496 0.61955
## DENSITYLow:SEASONWinter -0.8671 0.5338 -1.624 0.10432
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(9.0225) family taken to be 1)
##
## Null deviance: 183.269 on 41 degrees of freedom
## Residual deviance: 54.883 on 34 degrees of freedom
## AIC: 293.09
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 9.02
## Std. Err.: 3.69
##
## 2 x log-likelihood: -275.095
Conclusions:
tidy(quinn.nb, conf.int=TRUE)
Now if we exponentiate to put these estimates on the response scale:
tidy(quinn.nb, conf.int=TRUE, exponentiate=TRUE)
Conclusions:
If we compare the above to the overdispersed poisson, we see that the estimates are the same, but that the standard errors are underestimated and hence the confidence intervals are narrower (for overdispersed model)
tidy(quinn.glm, conf.int=TRUE, exponentiate=TRUE)
In order to tease appart the significant interaction(s), it might be useful to explore the effect of Density separately within each Season.
quinn.nb %>% emmeans(pairwise~DENSITY|SEASON)
## $emmeans
## SEASON = Spring:
## DENSITY emmean SE df asymp.LCL asymp.UCL
## High 2.303 0.187 Inf 1.935 2.67
## Low 2.416 0.200 Inf 2.024 2.81
##
## SEASON = Summer:
## DENSITY emmean SE df asymp.LCL asymp.UCL
## High 3.875 0.148 Inf 3.584 4.16
## Low 3.091 0.161 Inf 2.775 3.41
##
## SEASON = Autumn:
## DENSITY emmean SE df asymp.LCL asymp.UCL
## High 2.979 0.164 Inf 2.657 3.30
## Low 2.904 0.203 Inf 2.505 3.30
##
## SEASON = Winter:
## DENSITY emmean SE df asymp.LCL asymp.UCL
## High 1.735 0.219 Inf 1.306 2.16
## Low 0.981 0.402 Inf 0.192 1.77
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## SEASON = Spring:
## contrast estimate SE df z.ratio p.value
## High - Low -0.1133 0.274 Inf -0.413 0.6793
##
## SEASON = Summer:
## contrast estimate SE df z.ratio p.value
## High - Low 0.7836 0.219 Inf 3.577 0.0003
##
## SEASON = Autumn:
## contrast estimate SE df z.ratio p.value
## High - Low 0.0748 0.261 Inf 0.286 0.7749
##
## SEASON = Winter:
## contrast estimate SE df z.ratio p.value
## High - Low 0.7538 0.458 Inf 1.646 0.0999
##
## Results are given on the log (not the response) scale.
Conclusions:
Or we could express these on a response scale.
quinn.nb %>% emmeans(pairwise~DENSITY|SEASON, type='response')
## $emmeans
## SEASON = Spring:
## DENSITY response SE df asymp.LCL asymp.UCL
## High 10.00 1.87 Inf 6.93 14.44
## Low 11.20 2.24 Inf 7.57 16.58
##
## SEASON = Summer:
## DENSITY response SE df asymp.LCL asymp.UCL
## High 48.17 7.13 Inf 36.03 64.39
## Low 22.00 3.55 Inf 16.03 30.19
##
## SEASON = Autumn:
## DENSITY response SE df asymp.LCL asymp.UCL
## High 19.67 3.23 Inf 14.26 27.13
## Low 18.25 3.71 Inf 12.25 27.19
##
## SEASON = Winter:
## DENSITY response SE df asymp.LCL asymp.UCL
## High 5.67 1.24 Inf 3.69 8.70
## Low 2.67 1.07 Inf 1.21 5.87
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## SEASON = Spring:
## contrast ratio SE df z.ratio p.value
## High / Low 0.893 0.245 Inf -0.413 0.6793
##
## SEASON = Summer:
## contrast ratio SE df z.ratio p.value
## High / Low 2.189 0.480 Inf 3.577 0.0003
##
## SEASON = Autumn:
## contrast ratio SE df z.ratio p.value
## High / Low 1.078 0.282 Inf 0.286 0.7749
##
## SEASON = Winter:
## contrast ratio SE df z.ratio p.value
## High / Low 2.125 0.973 Inf 1.646 0.0999
##
## Tests are performed on the log scale
Conclusions:
quinn.nb %>% emmeans(pairwise~DENSITY|SEASON, type='response') %>%
confint
## $emmeans
## SEASON = Spring:
## DENSITY response SE df asymp.LCL asymp.UCL
## High 10.00 1.87 Inf 6.93 14.44
## Low 11.20 2.24 Inf 7.57 16.58
##
## SEASON = Summer:
## DENSITY response SE df asymp.LCL asymp.UCL
## High 48.17 7.13 Inf 36.03 64.39
## Low 22.00 3.55 Inf 16.03 30.19
##
## SEASON = Autumn:
## DENSITY response SE df asymp.LCL asymp.UCL
## High 19.67 3.23 Inf 14.26 27.13
## Low 18.25 3.71 Inf 12.25 27.19
##
## SEASON = Winter:
## DENSITY response SE df asymp.LCL asymp.UCL
## High 5.67 1.24 Inf 3.69 8.70
## Low 2.67 1.07 Inf 1.21 5.87
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## SEASON = Spring:
## contrast ratio SE df asymp.LCL asymp.UCL
## High / Low 0.893 0.245 Inf 0.522 1.53
##
## SEASON = Summer:
## contrast ratio SE df asymp.LCL asymp.UCL
## High / Low 2.189 0.480 Inf 1.425 3.36
##
## SEASON = Autumn:
## contrast ratio SE df asymp.LCL asymp.UCL
## High / Low 1.078 0.282 Inf 0.646 1.80
##
## SEASON = Winter:
## contrast ratio SE df asymp.LCL asymp.UCL
## High / Low 2.125 0.973 Inf 0.866 5.22
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
It might be useful to capture these pairwise contrasts so that we can graph them as a caterpillar plot.
In the following, I will present the effects on a log (based 2) axis. Such a scale presents doubling and halving (etc) equidistant from 1.
eff <- (quinn.nb %>%
emmeans(pairwise~DENSITY|SEASON, type='response') %>%
confint)$contrasts %>%
as.data.frame
eff %>%
ggplot(aes(y=ratio, x=SEASON)) +
geom_pointrange(aes(ymin=asymp.LCL, ymax=asymp.UCL)) +
geom_hline(yintercept=1, linetype='dashed') +
scale_x_discrete(name='') +
scale_y_continuous(name='Density effect (High vs Low)', trans=scales::log2_trans(),
breaks=scales::breaks_log(base=2)) +
coord_flip(ylim=c(0.25, 4)) +
theme_classic()
newdata = emmeans(quinn.nb, ~DENSITY|SEASON, type='response') %>% as.data.frame
head(newdata)
ggplot(newdata, aes(y=response, x=SEASON, fill=DENSITY)) +
geom_pointrange(aes(ymin=asymp.LCL, ymax=asymp.UCL, shape=DENSITY),
position=position_dodge(width=0.2)) +
theme_classic() +
theme(axis.title.x=element_blank(),
legend.position=c(0.01,1), legend.justification = c(0,1)) +
annotate(geom='text', x='Summer', y=70, label='*', size=7) +
scale_shape_manual(values=c(21, 22))
library(pscl)
quinn.zip <- zeroinfl(RECRUITS ~ DENSITY*SEASON | 1, data=quinn, dist='poisson')
#quinn.resid <- simulateResiduals(quinn.zip, plot=TRUE)
summary(quinn.zip)
##
## Call:
## zeroinfl(formula = RECRUITS ~ DENSITY * SEASON | 1, data = quinn, dist = "poisson")
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## -2.4510 -0.8645 0.1262 0.7225 2.8711
##
## Count model coefficients (poisson with log link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.3026 0.1291 17.836 < 2e-16 ***
## DENSITYLow 0.1133 0.1858 0.610 0.54191
## SEASONSummer 1.5721 0.1419 11.081 < 2e-16 ***
## SEASONAutumn 0.6763 0.1586 4.266 1.99e-05 ***
## SEASONWinter -0.5680 0.2147 -2.646 0.00815 **
## DENSITYLow:SEASONSummer -0.8970 0.2134 -4.202 2.64e-05 ***
## DENSITYLow:SEASONAutumn -0.1881 0.2381 -0.790 0.42957
## DENSITYLow:SEASONWinter 0.2165 0.4534 0.478 0.63300
##
## Zero-inflation model coefficients (binomial with logit link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.0037 0.7305 -4.112 3.92e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Number of iterations in BFGS optimization: 13
## Log-likelihood: -151.3 on 9 Df
#tidy(quinn.zip, conf.int=TRUE, exponentiate = TRUE)
exp(-3.0037)
## [1] 0.0496032
quinn.zip1 <- zeroinfl(RECRUITS ~ DENSITY*SEASON | SEASON, data=quinn, dist='poisson')
#quinn.resid <- simulateResiduals(quinn.zip)
summary(quinn.zip1)
##
## Call:
## zeroinfl(formula = RECRUITS ~ DENSITY * SEASON | SEASON, data = quinn,
## dist = "poisson")
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## -3.5327 -1.0591 0.1537 0.9216 3.6831
##
## Count model coefficients (poisson with log link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.3026 0.1291 17.836 < 2e-16 ***
## DENSITYLow 0.1133 0.1858 0.610 0.54191
## SEASONSummer 1.5721 0.1419 11.081 < 2e-16 ***
## SEASONAutumn 0.6763 0.1586 4.266 1.99e-05 ***
## SEASONWinter -0.5680 0.2147 -2.646 0.00815 **
## DENSITYLow:SEASONSummer -0.8970 0.2134 -4.202 2.64e-05 ***
## DENSITYLow:SEASONAutumn -0.1881 0.2381 -0.790 0.42958
## DENSITYLow:SEASONWinter 0.2291 0.4375 0.524 0.60045
##
## Zero-inflation model coefficients (binomial with logit link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.157e+01 1.454e+04 -0.001 0.999
## SEASONSummer -2.543e-08 2.013e+04 0.000 1.000
## SEASONAutumn -2.119e-08 2.107e+04 0.000 1.000
## SEASONWinter 2.031e+01 1.454e+04 0.001 0.999
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Number of iterations in BFGS optimization: 13
## Log-likelihood: -148 on 12 Df
exp(-3.0037)
## [1] 0.0496032
quinn.zinb <- zeroinfl(RECRUITS ~ DENSITY*SEASON | 1, data=quinn, dist='negbin')
AICc(quinn.zip, quinn.zinb)
summary(quinn.zinb)
##
## Call:
## zeroinfl(formula = RECRUITS ~ DENSITY * SEASON | 1, data = quinn, dist = "negbin")
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## -1.981e+00 -7.511e-01 -1.522e-05 5.289e-01 2.869e+00
##
## Count model coefficients (negbin with log link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.3026 0.1875 12.283 < 2e-16 ***
## DENSITYLow 0.1133 0.2742 0.413 0.67933
## SEASONSummer 1.5721 0.2389 6.580 4.69e-11 ***
## SEASONAutumn 0.6763 0.2492 2.714 0.00664 **
## SEASONWinter -0.5680 0.2881 -1.971 0.04869 *
## DENSITYLow:SEASONSummer -0.8969 0.3509 -2.556 0.01059 *
## DENSITYLow:SEASONAutumn -0.1881 0.3788 -0.496 0.61957
## DENSITYLow:SEASONWinter -0.8671 0.5338 -1.624 0.10432
## Log(theta) 2.1997 0.4091 5.378 7.55e-08 ***
##
## Zero-inflation model coefficients (binomial with logit link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -15.29 453.38 -0.034 0.973
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Theta = 9.0223
## Number of iterations in BFGS optimization: 43
## Log-likelihood: -137.5 on 10 Df
exp(-15.29)
## [1] 2.288956e-07
Quinn, G. P. 1988. “Ecology of the Intertidal Pulmonate Limpet Siphonaria Diemenensis Quoy et Gaimard. II Reproductive Patterns and Energetics.” Journalofexperimentalmarinebiologyandecology 117: 137–56.